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Abstract. We describe the spectral classification of white dwarfs and some of the physical 
processes important for their understanding. In the major part of this paper we discuss the 
input physics and computational methods for one of the most widely used stellar atmosphere 
codes for white dwarfs. 
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1 . Spectral classification of white 
dwarfs 

The classification scheme for white dwarfs 
(WD) developed in the beginning in analogy to 
the main sequence spectral types, with a distin- 
guishing letter "D" for degenerate object. DAs 
thus were stars with very strong Balmer lines, 
DBs had strong He I lines, DOs Hell. Today 
we know that this classification - in contrast to 
the main sequence - has not much to do with 
effective temperature, but is an indication of 
the photospheric composition. The classifica- 
tion used today was developed and described 
in detail in lSion et all (Il983i\ 

The main characteristic is the division into 
hydrogen-rich (DA) and helium-rich (DB, DO) 
atmospheres, but again in contrast to normal 
stars the most abundant element dominates 
with very few exceptions by several orders 
of magnitude. The explanation for this quasi 
mono-elemental composition is gravitational 
separation dSchatzmanl Il947i) . In the absence 
of significant competing macroscopic motions 
(stellar wind, meridional circulation, convec- 



Send offprint requests to: D. Koester 



tion) the heavier elements diffuse downward, 
leaving the lightest element present floating on 
top. The helium-dominated objects in this sce- 
nario must have lost their thin outer hydro- 
gen envelope during the formation phase of the 
white dwarf in the late stages of the asymptotic 
giant branch or planetary nebular phase. 

Besides the major types mentioned above 
one distinguishes DC (too cool to show 
any spectral feature, mostly helium-rich), DQ 
(atomic or molecular features of carbon), DAZ, 
DBZ, DZ (objects with traces of metals in 
hydrogen-rich or helium-rich atmospheres). 

The carbon in the DQ is assumed to be 
dredged up from deeper layers by the grow- 
ing conve ction zone in the superficia l he- 
lium l ayer dKoester et al.Lll982l:iPelletier et al.L 
Il986h . whereas the other heavy metals must be 
accreted from an outside source, either the in- 
terstellar matter, or some debris from a tidally 
disrupted asteroid. 

A spectral atlas showing many example 
spectra for a ll major ty pes has been published 
by .Wesemael et all ([T993i) . 
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2. Observational quantities 

Stellar parameters (effective temperature T^ff, 
surface gravity log g, abundances) are obtained 
from an analysis of spectroscopic or photomet- 
ric observations. If the surface of a star could 
be resolved, as for the sun, and if all relevant 
properties of our instrument were known, we 
could in principle determine the energy emit- 
ted by a small element of surface area, per unit 
time, wavelength interval, solid angle, into our 
line of sight. This quantity is called the (spe- 
cific) intensity, and in the case of a spherically 
symmetric star the only geometric variable for 
the surface value is the angle of emission rel- 
ative to the normal on the surface element &, 
that is 



/ = or / = I(fi) 



(1) 



with ju — cos &. If we cannot resolve the sur- 
face we can only measure the average intensity 
of the stellar disk I. More specifically, the en- 
ergy flux / arriving outside the terrestrial at- 
mosphere is related to this average intensity by 



(2) 



with average intensity 

7T/2 

7 = 2 J" /()?) cos()?) sin(i?) t/i? 




Un)ndn 



and the solid angle of the star 



Q 



D2 



(3) 



(4) 



with radius R and distance D. 

If we want to determine stellar parameters 
from a comparison of observed and theoreti- 
cally calculated spectra, the quantity which has 
to be calculated is thus the intensity / at the 
surface of the star The theory of stellar atmo- 
spheres has been developed by many authors 
over the past century and has reached a very 
mature state toda y. Clas sical works, still w orth 
reading, are e.g. lUnsold (.1968.) and .Mihalasl 



dl978h . ' 'Model atmospheres" and "synthetic 
spectra" as well as computer codes to calcu- 
late them are widely available. In the remain- 
der of this paper we will describe in detail the 
input physics and computational methods used 
by the author for his model atmospheres, which 
are used by many groups. 

The programming of the code was started 
by Dr. Thomas Gehren about 1975 with mi- 
nor contributions by myself. However, since 
then practically every routine has been com- 
pletely rewritten several times by the current 
author, and every remaining programming er- 
ror is only my fault. 

3. IVIodel atmosplieres and synthetic 
spectra 

The basic procedure is to specify the element 
abundances in the atmospheres, and the pa- 
rameters eff'ective temperature Teff and surface 
gravity log g, which are used as proxies for the 
"typical" values of the thermodynamic vari- 
ables temperature and pressure. Using a num- 
ber of simplifying assumptions and basic laws 
of physics this is sufficient to predict the radi- 
ation field (intensity) at the surface of the star. 
The most important assumptions are 

- homogenous, plane parallel layers: the 

depth of the atmosphere is considered to 
be very small compared to the radius of 
the star All matter quantities (density, pres- 
sure, temperature) depend only on one geo- 
metric variable, the height (in radial direc- 
tion) z- The intensity depends on z and the 
angle against the normal -t}, but not on the 
azimuthal angle. 

- hydrostatic equilibrium: at each point 
within the outer layers, which have a direct 
influence on the emerging radiation (i.e. the 
atmosphere or photosphere) the gradient of 
the gas pressure is in equilibrium with the 
gravitational attraction (plus possibly the 
transfer of momentum by photons). 

- radiative and convective equilibrium: 
there is no energy generation or loss within 
the atmosphere, only transport of the en- 
ergy generated in the deep interior. This 
transport can occur through radiation, heat 
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conduction, or convection; the total energy 
flux as determined by the parameter eff'ec- 
tive temperature is constant at all depths. 
- Local Thermodynamic Equilibrium: the 
matter is in thermal equilibrium corre- 
sponding to the local temperature at each 
layer, that is the ionization, excitation, dis- 
sociation of molecules etc. are governed 
by the usual relations of thermal equilib- 
rium (Boltzmann factors, Saha equation, 
Kirchhoff's law etc.). This is a very impor- 
tant assumption since it decreases the com- 
putational effort by several orders of 
magnitude. Please note that thermal equi- 
librium (i.e. the Planck function) is not as- 
sumed for the radiation field! Except for 
white dwarfs hotter than about 50000 K 
this LTE assumption is well justified. 

The code is divided in two major parts, 
which calculate in turn the physical structure 
of the outer layers (run of temperature, den- 
sity, pressure, absorption coefficients etc. with 
depth, this part will be called ATM here) and 
the surface intensities for many wavelengths 
(emerging spectrum, called SYN). There are 
auxiliary programs for additional necessary 
tasks, e.g. one for calculating the equation of 
state and absorption coefficients (KAPPA), an- 
other one for calculating equivalent widths of 
spectral lines or theoretical magnitudes in any 
photometric system (FILT), and so on. 

3.1. Equation of State (KAPPA) 

If no molecule formation has to be considered 
(e.g. high temperatures) and no elements be- 
sides H and He are present, the thermodynamic 
calculations are made directly in parallel with 
the determination of the atmospheric structure 
in the program ATM. Otherwise, these calcu- 
lations are made in KAPPA and the results 
(tables of matter density p, electron pressure 
Pe, entropy, absorption coefficients etc.) are 
stored in large two-dimensional tables as func- 
tion of temperature T and gas pressure Pg. The 
EOS, Saha equation for ionization, and dis- 
sociation equilibria for molecules are derived 
from a model Free Energy, which includes the 
ideal gas terms, Coulomb connections and an 



"Excluded Volume" term for the non-ideal in- 
teraction of neutral particles. Electron degener- 
acy is tested in all layers, but currently not im- 
plemented in the EOS, as it has been unimpor- 
tant in the range of parameters, where I have 
used my codes. 

Partition functions for H, Hel, and Hell 
are explicitly calculated using the lowest 
100 levels from die TOP BA SE database 
dCunto & Mendozal 119921: ICunto et al.L 
1 19931) and applying the occupation proba- 
bility w according t o the prescriptions of 
Humm er & MihalasI (|1988|); Mihalas et al] 
(il988h . For all other elements we use tables 
given bv lKuruc^ (11970 ) providing the partition 
function for a nominal cutoff 0. 1 e V below the 
ionization limit. The actual limit is calculated 
by a hydrogenic fit to the higher levels using a 
cutoff" determined from the non-ideal terms in 
the EOS. 

Currently dissociation equilibria are imple- 
mented for 20 molecules (H2, CH, NH, OH, 
MgH, SiH, CaH, C2, CN, CO, Nj, NO, O2, 
TiO, H7O HCN, HCO C^, CO?, N7O) us ing 
data from lKurucj(ll97 0') and Tatum (1966). 

The non-linear system of Saha and dissoci- 
ation equations together with the condition of 
neutrality and the definition of total gas pres- 
sure is solved with a Newton-Raphson itera- 
tion. 

3.2. Absorption coefficients (KAPPA) 

The absorption coefficient k describes the prob- 
ability w that a photon will interact (be ab- 
sorbed or scattered), when traveling a small 
distance ds in matter of density p 

w — Kpds — ^ n,- flj ds (5) 

For a dimensionless w, k has to have the dimen- 
sion of area per mass. It usually is the sum of 
many different interaction processes, with each 
contribution determined from the number den- 
sity of particles in the absorbing atomic state 
and the area a,, the cross section for this inter- 
action. The most important processes for white 
dwarfs and some sources of data or routines 
are (note that in most cases the data have been 
transformed by us and/or new routines written 
for our use): 
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bound-free and free-free absorption of 
neutral hydrogen: this can be calculated 
quasi-classically and corrected to quantum 
mechanics by the so-cal le d Gaunt factors 
dMenzel & PekerisL 1 19351: iKarzas & LatteR 
The free-free coefficient 



iKurucz 



for all other ions (with the exceptions noted 
below) is calculated hydrogen-like, 
bound-free and free-free absorption of the 

ion: numerical fits from John (1988.) . 
bound-free and free-free transitio ns of the 
H| ion : data and numerical fits from lBoggessI 

bound-free absorption of neutral helium: 

cross sections for the 43 lowest levels of Hel 
are taken from the TOPBASE database, 
free-free absorption of He : helium does not 
have a bound state as negative ion, so on ly the 
free-fr ee process is needed. Data are from lJohnI 
(11994 . 

bound-free and free-free absorption of the 
negative carbon ion C : data for the bf cross 
sect ions a re from Robinso n & GeltrnanI (Il967h 
and ICoop er & Martin ( 196^),^ for th e ff cross 
section from John & Williams ('1976'). 
bound-free transitions for elements other 
than H, He: these cross sections are mostly 
from the TOPBASE database if they are avail- 
able there, or hydrogen-like calculations other- 
wise. 

Thomson scattering by free electrons: the 

constant cross section per electron 



o- = 6.6527 X 10^25 cm" 



(6) 



is used. 

Rayleigh scattering by HI, Hel, H2 : cros s 
section fits are from Dalgamo (Il962l) : 
iDalgarno & W illiams ( 1962): lKuruc3(ll970l) . 
Molecular absorption: the calculations use 
the just-overlapping-line or smeared-line ap- 
proximation in the version developed by 
IZeidler-K T. & KoesteJ (|1982|) . This assumes 
that the density and broadening of rotational 
lines are so high that they form a quasi- 
continuum. Currently implemented are molec- 
ular data for C2, C3, and H2 molecules. 
Spectral line absorption: atomic data (ex- 
citation energies, oscillator strengths, line 
broadening constants) are obtained from a 
number of atomic databases, predominantly 



the Une lists from Kuruc z and collabora- 
tors (iKurucz & B ell'. '1995'), and the VALD 
(Vienna) database ( Kupka et a ll l2000i [T99^ 
Rvabchikova et al.L 119971: iPiskunov et al 



1995h . 



Satisfactory theories and data for the line 
profiles do exist f or the Stark b roadening of 
neutral hydrogen (iLemkel 1 19971: |Vidal et al^ 
I1973I) , and for 21 optical l i nes of neutral he- 
hum dBarnard et all 119691: iBeauchamp et al.L 
Il997h . These are so-called "unified theories", 
which attempt to describe the total line profile 
from core to the far wing. Similarly, the first 
three Lyman lines of H broadened by ionized 
and neutral perturbers and including a number 
of satellite features are well described by the 
work of Nicole A llard and collaborators (e.g. 
lAUardetaLl 12004 and many earlier papers). 

For all other processes the situation is 
much less satisfactory. Stark broadening pa- 
rameters for further Hel lines are pr ovided by 
iDimitriievic & Sahal-BrechotI d 19901) . In many 
later papers of the Belgrade group around 
Dimitrijevic similar data are provided for other 
elements. 

Below 8000 K for hydrogen-rich and 
16000 K for helium-rich atmospheres line 
broadening by neutral particles becomes im- 
portant. Since in most objects we have one 
dominating element, the interaction is usually 
between H-H or He-He. Resonance broaden- 
ing is thus important ( Ali & Griem, 1 1965b . as 
well as van der Waals interaction. On ly the 
first three Balmer lines (IBarklem et al.l I2OOO1) 



and some He transitions dLeo et al.Lll995h have 



so-called self-broadening theories; in the lat- 
ter case however for very low temperatures 
(300 K) only. These theories combine the ef- 
fects of resonance and van der Waals broaden- 
ing in a more consistent way. 

A few experimental measurements of 
broadening constants do exist, but in the vast 
majority of metal lines the Stark and van der 
Waals broadening constants can only be very 
roughl y estim ated by simple approximations 
(e.g. lUnsoldL Il968t ICowlevl Il97li: iGriemL 
Il966h . The line profile in the "impact approxi- 
mation" is then described by a Lorentz profile 
with these damping constants, which yields a 
Voigt profile after convolution with a Doppler 
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profile for the line-of-sight velocities of the 
emitting atoms. 



photons along the path ds, using the geometry 
described in Fig. [T| 



3.3. Atmospheric structure 

As mentioned above, the basic quantity for the 
description of the radiation field is the intensity 



I^I(z,A,ij) 



(7) 



with the geometrical height scale z measured 
from an arbitrary level outward of the star, 
wavelength A, and cosine of the angle against 
the z-axis ju = cos §. Useful quantities derived 
from this are the mean intensity (averaged over 
all directions, not to be confused with the disk- 
averaged intensity 7) 



(8) 



and the energy flux by radiation per unit area 
1 

F = 27rJ^ Ifid^ (9) 

-I 

Another useful quantity is 



(10) 



Note, that at the surface of a spherically sym- 
metric star with no radiation from the outside 
we have 



F ^nl 



(11) 



that is, the energy flux through the surface of 
a star is the quantity to be calculated for the 
comparison with non-resolved observations of 
a white dwarf. 



3.3.1 . The equation of radiative transfer 



dl 
dl 



psds — IpKds 



E ^ 

pKdz K 

dl . „ 
^57 



(12) 



e and k are the emission and absorption coeffi- 
cients per mass, t is a new depth variable called 
optical depth, replacing the geometric variable 
z, dT - -pKds. r = corresponds to the top of 
the atmosphere. S is the so-called source func- 
tion, S - eIk. 

For the solution of this differential equation 
two boundary conditions have to be specified. 
For the incoming radiation at the top 7(0, ^) - 
for < is usually assumed. At the bottom, 
at some large value » 1 the incoming in- 
tensity from below I(tb,ij) for fi > has to 
be specified. One possibility is to assume that 
at large depth the source function is the Planck 
function B (see below), expand it around the 
value tb and derive from the transfer equation 



I(T,fx)^B(T)+fl 



dB 
dr 



(13) 



or higher order approximations. 

For the absorption and emission coeffi- 
cients we have even at this phenomenologi- 
cal level of description to distinguish between 
two different processes. In the case of ab- 
sorption they are called "true absorption" 
and "scattering" cr, the corresponding emission 
processes are "thermal emission" e, and again 
scattering e,. Scattering means that a photon 
through interaction with matter changes its di- 
rection, but not the energy, whereas in the case 
of "true absorption" energy is absorbed by the 
matter and possibly re-emitted later with a dif- 
ferent energy; scattering processes thus do not 
lead to an energy coupling of the radiation field 
with matter. 

In our case of LTE Kirchhoff 's law states 



The equation of radiative transfer describes the e, 
balance between emission and absorption of 



Kt 



(14) 
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surface 




bottom 



Coordinate system for 
stellar atmospheres 
ds = dz / cos 9 



Fig. 1. Geometry for the radiative transfer equation 



an extremely powerful result. On the other 
hand, for isotropically distributed scattering 
particles we can derive 



e. - a- J 



(15) 



If both kinds of processes are important we can 
write the source function as 



^ _ e _ e, + 6, 
K K, + cr 



-B + 



cr 



Ki + cr 



Sometimes only the case with cr = and there- 
fore S = B is called LTE, or strict LTE. Since 
the inclusion of scattering in the source func- 
tion is not very difficult computationally, we do 
not make this distinction here. 

Because of the nature of the boundary con- 
ditions, given on both ends of the solution in- 
terval, the solution of the first-order equation is 
numerically difficult. Special precautions have 
to be taken to avoid exponentially increasing 
parasitic solutions. This is very ingenio usly 
avoided by the method of lFeautrie^ (Il964l) . We 
introduce new variables by dividing / into a 
symmetric (m) and an antisymmetric (v) part 
(we use fi > and write -/j. for the negative 
angles) 



m(t,//) = ^[/(t,ju) + /(t,-//)] 



v(T,yu) 



-[/(t,//)-/(t,-^)]. 



(17) 
(18) 



It is clear from the definitions that m(t, -fj) - 
u{T,fi) and v(t, -ju) - -v(T,/i), so we need the 



solution only for positive fi. Once u and v are 
known, we can always recover / as well as J 
and F. Writing the radiative transfer equation 
separately for positive and negative yU, adding 
and subtracting the two, we can derive the 
Feautrier equations 



-J (16) ^2l^ 



du 



(19) 
(20) 



The boundary conditions can easily be trans- 
formed to the new variables. The numerical 
solution of the second order transfer equation 
above is much easier and stable than for the 
first-order equation. 

3.3.2. Further constraints 

The solution of the transfer equation needs the 
values of B and absorption coefficients at each 
depth of the atmosphere, and therefore the ther- 
modynamic variables e.g. T and Pg. The ad- 
ditional constraints we have are the constant 
value of the transported energy flux and the hy- 
drostatic equation 



FUz) 



DO 



F{z, A) dA + F„,nv(z) = cr«r". (21) 
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Here cr^ is the radiation constant of Stefan's 
law and Famv the convective energy flux (see 
below). 

The balance between the gradient of the 
gas pressure, gravitational force and radiative 
force is 



dPg 1 r 

— =-pg+- Ki 

dz c J 



(A) F(A) dA 



(22) 



where the second term on the right side de- 
scribes the momentum transferred from the ra- 
diation field to the matter Defining a "stan- 
dard" absorption coefficient ks at a standard 
wavelength, or as a weighted mean over wave- 
length (e.g. the Rosseland mean), we can use 
the associated standard optical depth scale 
dxs - -pKsdz 

CO 

dP„ P If 

-r^^- K(A)FiA)dA. (23) 

dTs Ks CpKs J 



The hydrostatic equation thus provides a rela- 
tion between pressure scale, geometrical, and 
optical depth. For technical reasons we use the 
gas pressure as the independent variable, and 
the depths z and ts at each layer are derived 
quantities. 

Since p,Fconv, and the absorption coeffi- 
cients depend also on temperature, the typi- 
cal method of solution is to assume a temper- 
ature stratification T{Pg) and energy fluxes F 
(e.g. from a previous similar calculation or it- 
eration step) and solve the two constraint equa- 
tion above together with the radiative transfer. 
These equations together provide just enough 
equations for the unknowns, if the temperature 
structure is known. Since this is in general not 
the case, an iterative solution is necessary. The 
temperature dependent quantities are expanded 
around the current value, e.g. 

B{z, T, A) = B(z, To, A) + ^AT (24) 
dT 

The whole system of equations is then solved 
at once for the temperature corrections AT" and 
iterated with an improved temperature strat- 
ification, until the corrections become suffi- 
ciently small and all constraints are fulfilled. 



As can be seen from eqs.(21,22) the con- 
straints couple all wavelengths, which is re- 
sponsible for the huge number of unknowns. 
If the knowledge of the detailed angle depen- 
dence of / or M is not needed the computational 
burden in some intermediate steps can be con- 
siderably reduced by the m ethod of "variable 
Eddington factors" (lAuer & Mihalas. 1970). 
We start from the Feautrier equation eq.(19) 
and integrate over fi from to 1 



d^K 
'dr^ 



J-S 



(25) 



Under many conditions, in particular at large 
optical depths, the ratio K/J tends to a constant 
value 1 /3. We introduce a "variable Eddington 
factor" / - K/J to get 



dT^ 



^J-S. 



(26) 



Assuming / to be known, the structure of this 
equation is the same as that of the original 
Feautrier equation and can be solved with the 
same methods. The value of / has of course 
to be calculated from the original equation, but 
this can be done for one wavelength a time and 
therefore much faster. 

3.3.3. Convection 

Convection under the conditions of white 
dwarfs is highly turbulent. There is as yet 
no satisfactory theory describing the energy 
transport from first principles, nor any realistic 
numerical simulation, which could be imple- 
mented in routine calculations of atmospheric 
models. One has therefore to resort to the very 
crude mixing-length approximation, originally 
by Prandtl (1925 ), and adapted to stellar con- 
ditions by Bohm -Vitensd ([195 8). 

In the calculation of stellar evolution or 
even stellar atmospheres for "normal stars" our 
colleagues are generally content with one free 
parameter to describe the energy flux by con- 
vection in the mixing-length approximation. 
This parameter is the ratio of the mixing length 
to the pressure scale height a = I /Hp. In the 
case of white dwarfs we have gone further 
and use three numbers a,b,c, which appear in 
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the heuristic derivation of the theory, as ad- 
ditional free parameters. Different versions of 
the MLT are then denoted as e.g. MLl/a - 1, 
or ML2/a - 0.6, where MLl, ML2 describe 
the choice of a, b, c, and a the mixing - length 
jFbntaine et all 1 19811; iTassoul et all 1 19901: 
[Jordan et al.i. Il998h . Comparison of UV with 
optical spectra of ZZ Ceti DA white dwarfs 
around T^ff = 11000 - 12000 K has shown 
that a consistent descript ion is possible with 
MLl/2 (Hoe ster et alj 1 19941) or ML2/0.6 
dBergeron et al.Lll995i) . both of which describe 
what is called "intermediate efficiency" con- 
vection. The latter choice is at present used as 
quasi standard for DAs, also by this author. It 
is, however, quite clear that MLT in general is 
a very poor approximation and WD parameters 
are therefore still uncertain, when convection is 
important. 



3.3.4. Numerical solution 

The solution is obtained by a discretization of 
the depth scale Pg, the wavelengths A, and an- 
gles fi. Typical numbers for the grid points are 
4 values for /i between and 1, 100 depth 
points, and 1000 to 100000 wavelength points. 
Derivatives are approximated by difference 
quotients, and integrals by sums. For the inte- 
gration over angles (to obtain J,F) Gaussian 
quadratures are used for higher accuracy with 
few points. For the integration over wave- 
length simple trapezoidal rule or Simpson's 
rule are used. We then obtain a huge system 
of Unear equations for the variables / at each 
depth, wavelength, and angle, and the AT" at all 
depths. 

Fortunately the matrix of this system is 
very sparse - a tridiagonal band structure of 
sub-matrixes and some extra lir ies and columns 
from the constraint equations. iRybickil (Il97lh 
has demonstrated a very efficient elimination 
scheme, which results in a final linear system 
of rank equal to the number of depth points 
(e.g. typically 100), which is full and has to be 
solved by standard methods to determine the 
temperature corrections. When these correc- 
tions are deemed small enough (criteria used 
are often that the relative temperature correc- 
tions are smaller than 0.001 and the total flux 



at each depth is correct to 0. 1 percent), the at- 
mosphere structure is determined and all im- 
portant quantities (temperature, gas pressure, 
electron pressure, density, specific heat, adia- 
batic gradient, number densities of molecules, 
absorption coefficients) as function of depths 
are saved in a file for further use. 



3.4. Synthetic spectra (SYN) 

The calculation of the atmospheric structure 
with ATM needs of course also the radiation 
field, including the spectrum emerging from 
the surface. The reasons why we use a sepa- 
rate program SYN to calculate this again are 
the following: 

For the calculation of the atmospheric 
structure all wavelengths are coupled through 
the constraint equations, thus limiting the num- 
ber to typically a few 1000. On the other hand, 
because of the necessity to calculate the to- 
tal energy flux, the wavelength grid has to 
cover a large range from X-ray to far infrared. 
For the comparison with observations we typi- 
cally need only a smaller range, but with much 
higher wavelength resolution. As the structure 
is now known, we do not need to consider the 
wavelength coupling again, but can calculate 
the radiation field for each wavelength inde- 
pendently. 

Because of this reduced burden we are free 
to use many more (even weak) spectral lines, 
or more sophisticated line broadening theories. 
We can also include much more detailed calcu- 
lations of molecular absorption bands. 



3.4.1. Numerical method 

The emerging spectral energy distribution 
F(0, A) could be calculated using the Feautrier 
equations. However, for technical reasons we 
use a different method here. Integrating the 
original transfer equation over angle // we can 
derive an integral equation for the flux 



F(t) = 2n S(t')E2(t' -T)dT' 
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- 2;r J S(T')E2iT-T')dT'. (27) 



with the exponential integral function E2. In 
abbreviated form we write this as the flux in- 
tegral operator <!> 

F(T) = 0[5(r)]. (28) 

A special case is the flux emerging from the 
surface of the star, which as we know is equal 
to the observational quantity disk-averaged in- 
tensity F(0) = nl 



F{Qi) = 2n J SiT')E2iT')dT'. (29) 



In strict LTE S (r) = B{T{t)), which is known 
from the atmospheric structure, and the calcu- 
lation would be reduced to a simple integra- 
tion. In general, however, the source function 
may include a scattering term and we need the 
mean intensity J. Formally this can be derived 
directly from the transfer equation in a similar 
way as for F, with the result 

00 

J{T) = ^ (^') - ^'D (30) 


with the exponential integral E\. The integral 
operator in this equation is called the A opera- 
tor 



3.5. Theoretical photometry and 
equivalent widths 

The final results of SYN are stored in a binary 
disk file. This file contains a table of fluxes 
(or intensities) as a function of vacuum wave- 
lengths, the "synthetic spectrum" at the stellar 
surface, and in addition the basic parameters 
and structure data of the atmosphere model. 
These data are in a format, which can be used 
as input for the ATM program to start the iter- 
ation for a similar model. 

Auxiliary programs are available to trans- 
form these data, e.g. to air wavelengths, or 
into an ASCII file for use by other authors. 
These "export" files contain the flux table and 
a header similar to the FITS headers with 
all important parameters of the calculation. I 
strongly encourage my users to never separate 
this header from the table. 

A program FILT calculates equivalent 
widths of spectral lines from the flux table. It 
can also calculate theoretical photometry in ar- 
bitrary filter systems as e.g. 

CO 

y = -2.5 log J IiA)SviA)dA + Cv (32) 



with the total transmission 5 y of the filter plus 
optics, terrestrial atmosphere, etc. The constant 
Cv has to be determined from standard stars 
with known absolutely calibrated spectrum and 
measured magnitude in the corresponding sys- 
tem. Very often Vega is used for this purpose. 



/(r) = A[S (t)] = a [aB + (1 - a) J] . (31) 

We call this a formal solution, since S on the 
right hand side also contains the unknown J, 
which makes this an integral equation. 

For the numerical solution the depth scale 
is discretisized again, transforming the con- 
tinuous variables J,B into vectors and the 
A operator into a matrix. We use an 18- 
point Gaussian integration formula; the emerg- 
ing flux F(Q,A) = nl as well as the inten- 
sity 7(0, /l,ju) can finally be calculated from 
the source function 5 by a simple integration 
(summation). 



4. Some very technical remarks and 
outlook 

The code is currently written in the program- 
ming language FORTRAN??, but slowly - as 
time permits - transformed to FORTRAN95, 
which is much less prone to programming er- 
rors and much easier to maintain. Although 
considered a very old-fashioned programming 
language by many (who most likely never used 
it), 1 have been able to use my code over more 
than 30 yeiirs on dozens of computers and op- 
erating systems, in most cases without ever 
changing a single line of code. I am very grate- 
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ful for that and do not plan to ever change to 
another language. 

Since more than twenty years, and through 
very fundamental changes, the complete code 
has been under control of a version control 
system, starting with SCCS, RCS, and about 
two years ago changed to MERCURIAL. This 
means that I can recall the complete programs 
for any date or any version back to about 1985, 
IF I know the relevant data (version, compila- 
tion date, calculation date, etc.). These identi- 
fiers are written into every output file, includ- 
ing the ASCII format of the synthetic spectra, 
which are so widely distributed. 

These files also contain the information for 
some free parameters, like the mixing-length 
version used, or changes to the Hummer- 
Mihalas occupation probabilities. For this rea- 
son I urge users, to always keep the header with 
the spectrum table, such that the code version 
used can be identified in case of problems or 
questions. Unfortunately the system is not per- 
fect, since it does not keep track of a few data 
files, which have to be changed for different 
calculations, the most important being the file 
with the spectral line data. Different databases 
provide quite often different data for oscillator 
strengths or broadening constants. Depending 
on what I beheve at the time of calculation 
to be the most reliable values, these change 
from time to time, and 1 cannot always recon- 
struct what has been used after some years. 1 
am thinking how to solve this problem, but so 
far without result. 

From the programming aspect, as men- 
tioned above the code is slowly moved to 
FORTRAN95. It is already very modular, with 
many modules free of any side effects and be- 
ing reused unchanged in different programs. 
Programming the way I now know it should be 
was a hard learning experience over decades, 
and I am glad that FORTRAN95 supports al- 
most all my ideas and preferences much bet- 
ter than the older versions. The current aim is 
to make the whole program system very user- 
friendly to be able to put it in the public domain 
under a GPL or similar license in a few years. 
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